{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Comparing ParSpMatVec to Julia's sparse matvecs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Random Real Matrix: Test y = beta * y + alpha * A*x\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: Method definition getLaplacian(Int64) in module Main at In[54]:7 overwritten at In[56]:7.\n",
      "WARNING: replacing docs for 'getLaplacian :: Tuple{Int64}' in module 'Main'.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "getLaplacian"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "L,Lin,Lib = getLaplacian(n::Int)\n",
    "\n",
    "returns discrete 2D Laplacian, interior part, and boundary part, on regular mesh with n^2 cells.\n",
    "\"\"\"\n",
    "function getLaplacian(n::Int)\n",
    "    h  = 1/n\n",
    "    dx = spdiagm((-ones(n,1),ones(n,1)),0:1,n,n+1)/h\n",
    "    d2x = dx'*dx\n",
    "    L  = kron(speye(n+1),d2x) + kron(d2x,speye(n+1))\n",
    "   return L\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n=10000\t Base=0.0001\t ParSpMatVec=0.0001\t speedup=1.3283\n",
      "n=14125\t Base=0.0001\t ParSpMatVec=0.0001\t speedup=1.6935\n",
      "n=19953\t Base=0.0002\t ParSpMatVec=0.0001\t speedup=1.2694\n",
      "n=28184\t Base=0.0002\t ParSpMatVec=0.0003\t speedup=0.6741\n",
      "n=39811\t Base=0.0003\t ParSpMatVec=0.0002\t speedup=1.5173\n",
      "n=56234\t Base=0.0006\t ParSpMatVec=0.0003\t speedup=1.7774\n",
      "n=79433\t Base=0.0006\t ParSpMatVec=0.0005\t speedup=1.3409\n",
      "n=112202\t Base=0.0011\t ParSpMatVec=0.0007\t speedup=1.6192\n",
      "n=158489\t Base=0.0013\t ParSpMatVec=0.0010\t speedup=1.2669\n",
      "n=223872\t Base=0.0018\t ParSpMatVec=0.0017\t speedup=1.0566\n",
      "n=316228\t Base=0.0026\t ParSpMatVec=0.0020\t speedup=1.2822\n",
      "n=446684\t Base=0.0039\t ParSpMatVec=0.0032\t speedup=1.2020\n",
      "n=630957\t Base=0.0054\t ParSpMatVec=0.0049\t speedup=1.0955\n",
      "n=891251\t Base=0.0075\t ParSpMatVec=0.0065\t speedup=1.1540\n",
      "n=1258925\t Base=0.0126\t ParSpMatVec=0.0092\t speedup=1.3721\n",
      "n=1778279\t Base=0.0156\t ParSpMatVec=0.0149\t speedup=1.0464\n",
      "n=2511886\t Base=0.0235\t ParSpMatVec=0.0228\t speedup=1.0286\n",
      "n=3548134\t Base=0.0357\t ParSpMatVec=0.0300\t speedup=1.1888\n",
      "n=5011872\t Base=0.0604\t ParSpMatVec=0.0438\t speedup=1.3790\n",
      "n=7079458\t Base=0.0775\t ParSpMatVec=0.0628\t speedup=1.2336\n",
      "n=10000000\t Base=0.1044\t ParSpMatVec=0.0924\t speedup=1.1298\n"
     ]
    }
   ],
   "source": [
    "using ParSpMatVec\n",
    "using Base.Test\n",
    "\n",
    "# modify these parameters to test different cases\n",
    "testRandom    = false    # sprandn or 2D Laplacian\n",
    "testTranspose = true     # test A'*x or A*x\n",
    "nvec          = 1        # number of vectors to multiply \n",
    "numProcs      = 1        # number of processors (for ParSpMatVec)\n",
    "nrun          = 20       # number of runs\n",
    "\n",
    "if testTranspose\n",
    "    f1 = Base.LinAlg.At_mul_B!\n",
    "    f2 = ParSpMatVec.Ac_mul_B!\n",
    "else\n",
    "    f1 = Base.LinAlg.A_mul_B!\n",
    "    f2 = ParSpMatVec.A_mul_B!\n",
    "end\n",
    "\n",
    "if testRandom\n",
    "    getA = n-> sprandn(n,n,5/n)\n",
    "else\n",
    "    getA = n-> getLaplacian(round(Int64,sqrt(n)))\n",
    "end\n",
    "\n",
    "N     = logspace(4,7,21)\n",
    "time1 = zeros(length(N),nrun)\n",
    "time2 = zeros(length(N),nrun)\n",
    "\n",
    "for k=1:length(N)\n",
    "    n = round(Int64,N[k])\n",
    "    \n",
    "    A = getA(n)\n",
    "\n",
    "    x = rand(size(A,2),nvec);  x = x*10 - 5;\n",
    "    y = rand(size(A,2),nvec);  y = y*10 - 5;\n",
    "    \n",
    "    if nvec==1\n",
    "        x = vec(x)\n",
    "        y = vec(y)\n",
    "    end\n",
    "\n",
    "    alpha = 123.56\n",
    "    beta = 543.21\n",
    "    y2 = copy(y)\n",
    "    \n",
    "    # use julia's matvec\n",
    "    for j=1:nrun\n",
    "        tic(); \n",
    "            y2 = f1(alpha,A,x,beta,y2)\n",
    "        time1[k,j] = toq()\n",
    "    end\n",
    "\n",
    "    # use Fortran code\n",
    "    y3 = copy(y)\n",
    "    tic(); \n",
    "    for j=1:nrun\n",
    "        tic()\n",
    "        f2( alpha, A, x, beta, y3, numProcs ); \n",
    "        time2[k,j] = toq()\n",
    "    end\n",
    "    \n",
    "    # print results\n",
    "    @printf \"n=%d\\t Base=%1.4f\\t ParSpMatVec=%1.4f\\t speedup=%1.4f\\n\" n mean(time1[k,:]) mean(time2[k,:]) mean(time1[k,:])/mean(time2[k,:])\n",
    "    @test norm(y3-y2) / norm(y) < 1.e-12\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.5.0-dev",
   "language": "julia",
   "name": "julia-0.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.5.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
